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1.0  INTRODUCTION 

Analysis  of  aerodynamic  flows  is  used  to  support  the  development  and  evaluation  of 
weapons  systems  tested  in  the  ground  test  facilities  at  AEDC.  Aerodynamic  solutions  are 
used  in  experimental  test  planning,  test  execution,  and  post-test  analysis  to  increase  test 
facility  productivity,  upgrade  test  data  quality,  and  reduce  test  expense.  The  most  direct 
approach  to  the  computation  of  turbulent  viscous  flows,  which  are  the  types  of  flows  usually 
encountered  at  AEDC,  is  to  solve  the  Reynolds-averaged  Navier-Stokes  equations.  Critical 
flow  areas  such  as  separated  flows  and  shock/boundary  layer  interactions  may  be  computed 
in  this  manner,  thereby  aiding  in  the  evaluation  of  a  specific  test  article  or  facility. 

Algorithms  for  solving  the  Navier-Stokes  equations  are  extremely  complex  and  require 
extensive  computer  resources.  Therefore,  it  is  economically  expedient  to  explore  alternative 
solution  techniques  which  may  be  used  in  a  flow  evaluation  without  large  computer  resource 
requirements.  An  inexpensive  approach  to  the  problem  is  by  means  of  a  component  method 
based  on  the  assumption  that  the  boundary-layer  equations  adequately  describe  the  viscous 
layer  near  a  wall  or  model,  and  that  the  exterior  flow  may  be  described  by  the  inviscid 
equations.  These  two  regions  of  flow  may  then  be  matched  at  a  common  boundary,  usually 
the  boundary- layer  edge. 

This  component  type  of  flow  solver  has  been  shown  to  be  an  effective  tool  in  the 
computation  of  viscous-inviscid  interacting  flows  as  cited  by  Lock  (Ref.  1),  Melnik  (Ref.  2), 
and  Le  Balleur  (Ref.  3).  This  method  was  demonstrated  for  two-dimensional  and 
axisymmetric  bodies  with  small  separated  flow  regions  in  Ref.  4.  Analysis  of  separated  flows 
may  be  accomplished  by  computing  the  boundary-layer  equations  using  an  inverse  method. 
A  review  of  computational  methods  capable  of  separated  flow  solutions  is  given  by  Le 
Balleur  (Ref.  3).  An  extension  of  inverse  component  methods  to  three-dimensional  problems 
is  feasible  and  has  the  potential  for  computing  separated  viscous  flows  over  complex  bodies 
with  significantly  less  computer  resources  than  Navier-Stokes  methods. 

The  component  method  may  be  referred  to  as  the  viscous-inviscid  interaction  procedure. 
Calculation  of  the  inviscid  exterior  flow  field  can  be  obtained  by  solution  of  the  Euler 
equations  which  approximate  the  Navier-Stokes  equations  as  the  Reynolds  number 
approaches  infinity  (and  removes  all  energy  dissipation  processes  from  the  flow).  The 
interior  region  may  be  computed  with  the  boundary-layer  equations  obtained  by  imposing 
the  Prandtl  limit  to  the  Navier-Stokes  equations  and  may  be  approximately  solved  by  the 
integral  method.  The  exterior  and  interior  regions  about  an  airfoil  are  depicted  in  Fig.  1. 

A  steady  interaction  code1  (one  in  which  the  boundary-layer  solution  method  is  not  time- 
dependent)  exists.  It  employs  a  two-dimensional  integral  compressible  boundary- layer 


1  The  terms  “steady”  and  “unsteady”  interactions  refer  to  the  boundary  layer  portion  of  the  code  only.  The 
Euler  code  is  time-dependent  for  both  cases. 
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Figure  1.  Flow  regions  about  an  airfoil. 


method,  developed  by  Whitfield  (Ref.  5),  in  conjunction  with  a  three-dimensional,  time- 
dependent  Euler  code  (Ref.  6).  The  two-dimensional  boundary-layer  method  can  be  applied 
to  a  three-dimensional  wing  by  using  strip  theory  when  the  crossflow  is  small  (since 
boundary  layers  are  sensitive  to  small  crossflow  pressure  gradients).  Calculations  of  two- 
dimensional  boundary  layers  are  made  at  various  spanwise  cross  sections.  However,  there 
exists  a  need  for  a  fully  three-dimensional  viscous  interaction  method  to  be  used  for  finite 
wing  applications.  This  would  allow  flows  to  be  analyzed  with  the  boundary  layers 
computed  over  complete  wings. 

There  are  problems  associated  with  a  three-dimensional  viscous-inviscid  interaction 
code.  In  two-dimensional  boundary-layer  codes,  the  solution  process  marches  in  the  flow 
direction  and  computes  boundary-layer  quantities  at  each  spatial  station.  In  three- 
dimensional  problems  the  solution  process  is  restricted  to  one  direction  but  there  are  two 
computational  directions  (e.g.,  the  chordwise  and  spanwise  directions  on  a  finite  wing).  If 
this  problem  is  alleviated  by  aligning  the  computational  mesh  on  the  body  with  the  flow 
direction,  the  mesh  becomes  flow-dependent  and  must  be  altered  after  each  update  to  the 
flow.  This  makes  any  solution  process  both  difficult  and  lengthy  because  of  the  constant 
altering  of  the  mesh  between  the  inviscid  code  and  the  boundary-layer  code.  Therefore,  it 
would  be  useful  to  produce  a  boundary-layer  method  which  utilizes  the  spacing  of  the 
inviscid  mesh  and  still  provides  a  solution  on  a  three-dimensional  body.  This  can  be 
accomplished  if  the  boundary-layer  equations  include  the  time-dependent  terms  and  the 
entire  boundary-layer  on  the  body  is  solved  at  each  instant  in  time.  The  solution  process 
could  march  by  means  of  a  time-stepping  parameter  until  a  steady-state  solution  is  obtained. 

However,  before  three-dimensional  flow  solutions  are  attempted  by  the  unsteady 
boundary- layer  method,  attention  should  be  focused  on  simpler  two-dimensional  problems 
solved  by  the  same  method.  This  allows  exploration  of  the  characteristics  of  the  unsteady 
formulation  without  the  complications  associated  with  another  added  dimension.  The 
unsteady  results  can  be  validated  by  comparisons  with  solutions  obtained  from  the  two- 
dimensional  steady  interaction  code. 

The  intention  of  this  work  is  to  produce  the  unsteady  interaction  code  and  compare  its 
solution  with  that  of  the  steady  interaction  code  and  with  data. 
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2.0.  DERIVATION  OF  THE  UNSTEADY  INTEGRAL  BOUNDARY -LAYER 
EQUATIONS  FOR  TWO-DIMENSIONAL  COMPRESSIBLE  FLOW 

The  inner  region  solution  is  constructed  from  the  integral  form  of  the  boundary-layer 
equations.  These  equations  are  obtained  via  the  Reynolds-averaged  Navier-Stokes  equations 
and  the  Prandtl  limit.  Thus,  when  expressed  in  conservation  form,  the  Navier-Stokes 
equations  are  (Ref.  7): 


where 


and 


are  the  dependent  variable  matrices  and 


Tij  =  £ 


( 


dvj 

3xj 


+ 


dvk 

X - - 

3xk 


k 


(1) 


(2) 


(3) 


(4) 


is  the  viscous  stress  tensor  expressed  in  terms  of  the  Lame  constants  p  and  X  (usually  referred 
to  as  the  first  and  second  coefficients  of  viscosity,  respectively)  and  is  the  Kronecker  delta. 
These  constants  are  related  to  the  bulk  coefficient  of  viscosity,  ?/,  by  the  equality,  ij  =  X  + 
2/i/ 3,  and  this  quantity  is  placed  equal  to  zero  by  adopting  Stokes’  hypothesis.  The  equation 
of  state  is  taken  to  be 


e  =  -E-  +  4  e  viVi  (5) 

7-1  2 

Finally,  q,  =  -kST/dx,  is  the  heat  flux  component  in  the  X;  direction  and  k  is  the  thermal 
conductivity.  The  form  of  Eq.  (1)  remains  intact  under  the  Euler  limit  operator  (E  such  that 
fi  —  0,X  —  0,  k  —  0),  but 
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E(fij)  - 


(6) 


is  the  limit  form  of  fjj,  and  all  dissipation-causing  terms  vanish. 

Before  the  Euler  limit  operator  is  applied  to  Eq.  (1)  it  is  necessary  to  introduce  the 
Reynolds  decomposition  in  the  form: 


</>=$+ 


(7) 


where 


<t>  =  E  (<«  and  E  (*')  =  0,  (8) 

with  E  an  averaging  operator.  The  Reynolds-averaged  equations  then  follow  when  the 
operator  E  is  applied  to  the  decomposed  equations.  All  nonlinear  terms  give  rise  to 
correlations  of  variable  fluctuations.  For  example,  a  velocity  product  —  which  arises  in  the 

f  t 

inertia  term  of  the  momentum  equation  —  will  produce  <  pv^vj  >  correlations  whose 
dimensional  form  allows  them  to  be  referred  to  as  Reynolds  stresses. 

The  Euler  equations  are  recovered  (as  the  appropriate  outer  flow  equations)  when  the 
Euler  limit  operator  E  is  applied  to  the  Reynolds-averaged  equations  if  it  can  be  assumed 
that  all  correlation  terms  vanish  in  the  far  field.  This  need  not  be  true  for  all  flows  of 
interest:  examples  to  the  contrary  include  wind  tunnel  flows  and  atmospheric  flows.  It  is  well 
known  that  free-stream  turbulence  can  have  a  significant  effect  on  the  development  of  a 
turbulent  boundary  layer  and  this  fact  may  well  impair  the  comparison  between  wind  tunnel 
data  and  theoretical  calculations  based  upon  the  classical  Euler  equations  for  the  external 
flow.  Nevertheless,  the  present  study  neglects  the  correlation  terms  in  the  Euler  limit 
equations. 

The  inner  region  equations  are  obtained  from  the  Reynolds-averaged  equations  by  use  of 
the  Prandtl  limit  operator  P,  where  P  is  a  product  operator  such  that: 

1.  The  coordinate  (x2,  say)  normal  to  the  wall  is  scaled  as  £  =  x2/  \fv  for  laminar  flows. 
Such  scaling  is  usually  adopted  for  turbulent  flow  without  justification. 

2.  The  operator  E  is  applied. 

The  result  of  this  combined  operation  is  the  pair  of  equations  (now  expressed  in  planar 
coordinates,  where  u  and  v  are  the  velocity  components  in  the  x-  and  y-directions, 
respectively): 
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3q_ 

dt 


+ 


d_ 

dx 


(eu)  + 


—  (ev)  -  o 

ay 


Q 


du 

at" 


+  eu 


au 

dx 


+  ev 


du 

ay 


dp  dr 

~r—  +  — —  » 
dx  dy 


(9) 

(10) 


in  which 


T 


<eu'v'  > 


01) 


represents  the  total  shear  stress  component  which  has  the  major  impact  on  a  two 
dimensional  boundary  layer. 


These  partial  differential  equations  may  be  transformed  into  ordinary  differential 
equations  in  space  by  multiplying  Eq.  (9)  by  um  +  V(rn+l),  and  Eq.  (10)  by  u"\  and 
summing  the  two  equations  to  form  the  mth-moment  equation: 


where  m  =  0  generates  the  momentum  equation,  and  m  =  1  creates  the  mean  flow  kinetic 
energy  equation.  The  momentum  equation  valid  at  the  edge  of  the  boundary  layer  has  been 
incorporated  to  replace  the  pressure  gradient,  dp/dx,  but  does  not  imply  the  assumption  of 
isentropic  flow. 


The  integral  momentum  equation  is  obtained  by  integrating  Eq.  (12)  from  0 
setting  m  =  0  and  introducing  the  integral  parameters 


F(l  — ) dy, 

o  \  0euc  / 


£  y<oo. 


(13) 


F-^fl  --V 

o  6eue  \  ue  } 


(14) 
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(15) 


To 

V2  GeUe 


06) 


in  which  5*  is  the  displacement  thickness  whose  significance  will  be  discussed  in  Section  5.0; 
9  is  the  momentum  thickness  or  the  distance  that  the  momentum  defect  extends  from  the 
body;  0eis  the  density  thickness;  and  Cf  is  the  local  skin  friction  coefficient  with  t0  = 
fidu/dy  |  y=0  since  the  Reynolds  stress  vanishes  on  the  wall.  Then  by  applying  the  boundary 
conditions 


u  =  0;v  =  0aty  =  0 


u  —  ue;  7  —  0  as  y  —  00 


the  integral  momentum  equation  becomes 


© 


1  a 

GeU« 


(e«u^) 


ue  dx 


(17) 


The  source  term  (A)  derives  from  the  time-dependent  terms  in  Eqs.  (9)  and  (10)  and  has 
two  components  involving  b*  and  0Q  which  come  from  both  the  continuity  and  momentum 
equations.  Term  B  describes  the  rate  of  change  of  momentum  thickness  flas  the  boundary 
layer  develops  under  the  imposed  pressure  gradient.  The  pressure  gradient  is  contained  in 
term  (c)  (here  expressed  as  an  external  velocity  gradient  where,  as  pointed  out  above,  the 
flow  has  not  been  assumed  to  be  isentropic).  Term  (c)  thus  represents  the  external  flow 
boundary  condition.  Finally,  term(B)  represents  the  only  contribution  of  the  shear  stress  t 
to  the  momentum  integral  equation  and  contains  only  the  boundary  value  tq  =  pdu/dy  |y=0. 
The  skin  friction  coefficient  cf  is  included  empirically  and  is  usually  expressed  in  the  form  cf 
=  c^Refl,  H,M),  as  discussed  in  Section  3.0. 


The  mean  flow  kinetic  energy  equation  is  obtained  by  setting  m  =  1  in  Eq.  (12), 
integrating  from  0  <  y  <  00  and  introducing  the  integral  parameters 
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(18) 


(19) 


(20) 


in  which  5*  is  the  velocity  thickness;  8*  is  the  energy  dissipation  thickness;  and  Dg  is  the 
shear  work  integral  term.  Then  by  applying  the  previous  boundary  conditions,  the  mean 
flow  kinetic  energy  equation  (or  first  moment  equation)  becomes 


1 _ 

2&Ue  dt 


eeu  l(6  +  8*-0ll) 


+  ^0 

Ue  3t 


+ 


1 

2eeue3 


d_ 

dx 


(Qeufr) 


|  3*-3u  3ue  _  CfDg  _  q 


Ue 


dx 


(21) 


©  ® 


Term  (^)  is  again  a  source  term  and  term  (c)  involves  the  external  pressure  gradient. 
Term  (5)  is  the  rate  of  change  of  mean  flow  kinetic  energy  thickness  8*.  Term  (3), 
containing  the  shear  work  integral  De,  is  the  most  interesting  term  since  this  is  the  major 
contribution  from  the  Reynolds  shear  stress  <pu'v'>.  From  Eq.  (20)  it  is  seen  that  the 
shear  work  integral  involves  an  integral  of  the  quantity  rdu/dy  across  the  boundary  layer. 


If  Eqs.  (17)  and  (21)  are  formed  for  the  three-dimensional  problem,  several  differences 
would  be  evident.  Two  momentum  equations,  each  with  three  velocity  terms,  would  be 
required  in  addition  to  the  three-dimensional  unsteady  continuity  equation.  Also,  a  second 
shear  stress  term  would  be  required  for  the  spanwise  direction  on  a  wing.  Two  sets  of 
integral  parameters  would  also  be  required;  one  in  each  direction  of  the  airfoil  —  chordwise 
and  spanwise. 

Finally,  it  can  be  noted  that  Eqs.  (17)  and  (21)  contain  seven  unknowns  15*,0,0*,0o,6*,Cr, 
DEj  and  can  only  be  solved  when  five  additional  equations  are  available.  Experimental 
correlations  are  used  in  this  capacity  in  the  present  work,  and  will  be  discussed  in  Section 
3.0. 
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3.0  STEADY  AND  UNSTEADY  BOUNDARY -LAYER  METHODS 
3.1  STEADY  BOUNDARY-LAYER  METHODS 


Many  of  the  basic  ideas  used  in  the  steady  formulations  are  maintained  in  the  unsteady 
case.  A  steady  compressible  turbulent  integral  boundary-layer  code  was  developed  by 
Whitfield,  Ref.  5,  to  produce  quick  and  accurate  solutions  to  two-dimensional  boundary 
layers  on  adiabatic  walls.  The  input  required  is  a  distribution  of  pressure  on  the  body  which 
is  usually  obtained  from  experimental  measurements  in  the  present  applications.  (This 
method  is  referred  to  as  the  direct  method.  Conversely,  a  method  using  a  5*  distribution  as 
input  is  termed  an  inverse  method.)  Whitfield’s  method  solves  the  integral  x-momentum  and 
mean  flow  kinetic  energy  equations  formulated  from  the  differential  continuity  and 
x-momentum  equations.  The  integral  equations  solved  are  (from  Ref.  5). 


1 

QeUc 


-J-  (QcU&) 
dx 


5*  duE 
+ - - 

ue  dx 


if 
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(22) 


1 

2ecUe 


-J-  (eeUe0*) 
dx 
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5* -a;  due 
uE  dx 


(23) 


These  equations  differ  from  the  unsteady  equations  developed  in  Section  2.0  only  by  the 
time  derivative  terms  [the  A  terms  from  Eqs.  (17)  and  (21)].  The  number  of  unknowns  in 
Eqs.  (22)  and  (23)  is  reduced  to  six  since  0C  only  occurs  in  the  unsteady  terms.  The 
correlations  obtained  later  will  be  used  as  the  remaining  equations  required  to  solve  the 
system, 


Whitfield’s  approach  to  solution  of  Eqs.  (22)  and  (23)  was  to  establish  various  shape 
factors: 


Ha*  =  &*/6,  .  (24) 

He.  =  0*/6,  (25) 

H5‘  =  K/0,  (26) 

and  correlate  them  with  H  and  the  boundary-layer  edge  Mach  number  by  numerically 
integrating  (Simpson’s  rule)  the  velocity  profile  of  Whitfield  and  the  expression  for  q/qc 
(Ref.  5).  The  correlations  relate  compressible  and  incompressible  flows  by  means  of  the 
Mach  number  dependence  included  in  the  correlations.  In  addition,  the  skin  friction  (Cf)  was 
correlated  with  H  and  Re*  by  using  the  work  of  Winter  and  Gaudet  (Ref.  8),  who  established 
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the  relationship  between  low-speed  skin  friction  (Cf)  and  the  local  skin  friction  (Cf)  in  a 
compressible  flow.  Their  relation  is 


■V  =  (1  +  JrT  (27) 

in  which  Me  is  the  Mach  number  at  the  boundary-layer  edge.  Initially,  this  relationship  was 
established  for  adiabatic  walls,  zero  pressure  gradient,  and  air  as  the  working  fluid,  but  it 
was  hypothesized  (Ref.  5)  that  the  relationship  is  also  true  if  the  appropriate  expression  for 
Cf  is  used  in  instances  in  which  the  pressure  gradient  does  not  vanish.  A  relationship  for  cf 
was  extracted  from  White  {Ref.  9,  Eq.  (6.180b)]  where  Cf  is  a  function  of  H  and  Re#;  that  is, 

cf  =  0.3  e<  - 1-33»>  (log  Re#)  <  - 1 -74-0.31  hi  (28) 

The  correlation,  Eq.  (28),  assumes  a  smooth  wall.  Effects  of  roughness  could  be  included  by 
introducing  new  parameters  into  the  correlation. 

The  shear  work  integral  term,  DE,  is 

De  =  ? 

O 

which  is  basically  a  turbulent  energy  production  term  which  describes  the  conversion  of 
mean  flow  energy  into  turbulent  energy.  Equation  (29)  was  obtained  by  direct  integration 
employing  the  Cebeci-Smith  two-layer  eddy  viscosity  turbulence  model  (Ref.  10)  for  r  and 
Whitfield’s  turbulent  velocity  profile  for  du/9y.  The  Cebeci-Smith  turbulence  model  was 
chosen  primarily  because  of  its  simplicity  and  computational  speed.  However,  the  shear 
work  integral  term  was  later  numerically  correlated  with  H,  edge  Mach  number,  Me,  and  Re# 
as  reported  in  Ref.  11.  As  a  result,  all  terms  in  the  steady  integral  equations  [Eqs.  (22)  and 
(23))  are  either  input,  or  may  be  determined  by  the  correlations.  The  integral  equations  are 
then  solved  by  the  predictor-corrector  method  of  Nash  (Ref.  12).  Results  from  Whitfield’s 
code  (Ref.  5)  show  acceptable  agreement  with  a  multitude  of  boundary-layer  data.  The 
program  is  known  by  the  acronym  SWIM  (Shear  Work  Integral  Method). 

The  SWIM  code  provides  a  reasonable  prediction  of  the  flow  separation,  if  separation  is 
designated  as  the  condition  at  which  the  numerical  algorithm  breaks  down.  The  problem 
arises  primarily  because  the  boundary-layer  equations  are  parabolic  in  nature  and  exhibit  a 
singularity  at  separation  (Ref.  13).  Physically,  this  means  the  influence  of  the  flow  is  felt 
only  in  the  downstream  direction.  Therefore,  after  separation  occurs,  the  flow  information 
cannot  influence  the  flow  upstream,  producing  a  singularity  at  the  separation  point. 
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Because  of  the  importance  of  separated  flows,  work  was  needed  to  extend  Whitfield’s 
method  to  the  separated  flow  regime.  First,  improvements  to  Whitfield’s  velocity  profile 
correlation  to  include  reverse  flows  were  required.  Swafford  (Ref.  14)  made  improvements 
and  succeeded  in  extending  Whitfield’s  shape  factor  correlations  and  the  skin  friction 
correlation  of  White  to  the  separated  flow  region. 

Second,  the  singularity  at  the  separation  point  needed  to  be  eliminated.  Catherall  and 
Mangier  report  (Ref.  15)  that  the  singularity  may  be  removed  by  the  specification  of  the 
displacement  thickness  (5*)  or  the  wall  shear  distribution  rather  than  a  pressure  distribution 
as  in  the  SWIM  code.  The  displacement  thickness  may  be  prescribed  when  the  boundary 
layer  is  computed  by  the  inverse  method.  Swafford’s  correlations  reveal  that  when  H 
becomes  large  enough  (—3.2),  a  separated  flow  profile  and  a  negative  skin  friction  will 
emerge.  A  plot  of  the  skin  friction  correlation  is  shown  in  Fig.  2.  Because  of  the  limited 
amount  of  experimental  data  available  for  separated  flows,  the  correlations  used  in  the 
separated  regions  cannot  yet  be  considered  universally  applicable,  especially  for  separated 
flows  in  which  H  >  4.  However,  the  correlations  are  generally  considered  good  for  H  <  4. 
The  Cebeci-Smith  turbulence  model  is  retained  for  the  separated  region  with  the  assumption 
that  the  turbulence  model  is  valid  in  reverse  flow.  This  is  questionable  since  the  turbulent 
energy  production  term  is  very  small  in  a  separated  flow  region.  The  outer  portion  of  the 
boundary  layer  is  still  “wake-like,”  however. 

The  SWIM  code  has  been  altered  (Ref.  4)  to  make  use  of  the  inverse  method  and  was 
coupled  with  an  Euler  equation  code  to  produce  a  steady,  inverse  viscous-inviscid 
interaction  code.  The  interaction  code  works  well  for  both  attached  and  separated  flow  and 
many  two-dimensional  problems  may  now  be  solved.  The  interaction  process  will  be 
discussed  in  more  detail  in  Section  5.0. 

3.2  UNSTEADY  BOUNDARY-LAYER  METHODS 

An  unsteady  boundary-layer  method  based  upon  the  principles  used  in  the  SWIM  code 
and  the  separated  flow  correlations  was  devised  by  Whitfield.  It  solves  the  two-dimensional 
integral  compressible,  unsteady  boundary-layer  equations  derived  in  Section  2.0.  Time  is  an 
incremental  stepping  parameter  in  the  formulation.  Thus,  the  entire  boundary  layer  is 
determined  after  each  time  step.  It  must  be  emphasized  that  the  time-advancing  routine  is 
only  a  means  of  arriving  at  a  final  steady-state  solution.  Any  solution  in  the  interim  is  not 
representative  of  the  real  time  flow.  Hence,  the  temporal  gradients  are  not  necessarily 
correct. 

The  unsteady  integral  boundary-layer  equations  were  formulated  to  be  solved  by  the 
direct  method  (see  Appendix  A),  in  which  surface  pressures  are  used  as  input  to  the  code. 
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H 

Figure  2.  Correlation  of  incompressible  skin  friction,  Cf,  for 
attached  and  separated  flows. 

The  parameters  computed  from  the  solution  process  are  the  momentum  thickness  (0)  and  the 
incompressible  shape  factor  (H).  A  flat  plate  approximation  of  the  B  and  H  distributions  is 
used  to  begin  the  solution  process.  An  additional  shape  factor  (H#p)  correlation  was  required 
to  solve  the  unsteady  problem  and  was  obtained  in  the  same  manner  as  the  previous  shape 
factor  correlations.  After  each  time  step  of  the  unsteady  method,  a  new  distribution  of  9  and 
H  is  determined  on  the  entire  body.  The  new  distributions  are  used  as  an  approximation  for 
the  next  time  step.  The  iteration  process  is  continued  until  steady  state  is  obtained.  All  other 
integral  parameters  may  be  computed  by  means  of  the  shape  factor  correlations  once 
distributions  of  Mach  number  and  H  are  known.  The  Mach  number  is  known  from  the 
input.  The  unsteady  direct  solution  process  is  shown  in  Table  ]. 

Cousteix,  et  al.  (Ref.  16)  have  suggested  that  singularities  are  possible  with  an  unsteady 
direct  method.  Whitfield  developed  an  inverse  unsteady  method  to  avoid  such  singularities, 
as  was  done  in  the  steady  version.  The  unsteady  boundary-layer  equations  were  formulated 
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Table  1.  Solution  Procedure  for  the  Unsteady,  Direct 
Boundary-Layer  Method 

Step  Procedure 

L.  Input  initial  values  of  H(x)  and  0(x)  distributions  (flat  plate 
approximations)  and  Me(x). 

2.  Evaluate  all  other  parameters  using  H,  Me,  and  Re*  at  every  spatial 
station. 

3.  Solve  equations  for  88/ dt  and  9H/3t  at  every  spatial  station. 

4.  Compute  8,  H  at  next  time  step  (backward  difference  in  time)  at  all 
spatial  stations. 

5.  Use  8,  H  distributions  for  next  time  step. 

6.  Go  to  step  2  and  do  until  convergence  (where  86/ dt  and  dH/ dt  —  0). 

for  the  inverse  method  (see  Appendix  B)  in  which  a  8*  distribution  is  used  as  input.  The 
parameters  to  be  determined  by  the  time-dependent  inverse  solution  process  are  8  and  the 
edge  velocity  (ue).  As  was  done  in  the  direct  unsteady  method,  the  two  boundary-layer 
parameters  (8  and  ue)  which  are  obtained  in  the  solution  are  updated  in  time  until  steady 
state  is  reached.  The  unsteady  inverse  process  is  shown  in  Table  2. 

Whenever  a  5*  distribution  is  not  available,  an  inverse  boundary-layer  code  may  be  used 
in  a  pseudodirect  mode.  Here,  a  pressure  distribution  is  used  as  input  and,  by  the  iterative 
process  of  Carter  discussed  in  Section  3.0,  an  initial  6*  distribution  (using  flat  plate  data,  for 
example)  may  be  updated  until  convergence.  With  each  6*  distribution  during  the  solution 
process,  a  pressure  distribution  is  obtained  from  the  boundary-layer  code  and  compared  to 
the  input  value  of  pressure.  The  ratio  of  computed  pressures  to  input  pressures  is  used  in 
Carter’s  method  to  update  the  next  5*  distribution.  In  this  manner  a  pressure  distribution  is 
used  as  input,  yet  the  inverse  method  is  used  to  obtain  a  solution.  This  method  is  valid  only 
for  attached  flows.  It  is  not  known  why  this  method  will  not  compute  separated  flows  and 
an  inverse  method  will.  The  pseudodirect  method  was  used  in  the  present  work  when 
comparisons  between  direct  and  inverse  codes  were  made  for  attached  flows. 

Both  the  direct  and  inverse  methods  employ  a  Fourth-order,  variable  siep  Runge-Kutta 
scheme  in  the  solution  processes.  This  method  had  the  best  convergence  record  with  an 
upstream  spatial  differencing  scheme.  Central  differencing  was  tried  earlier  and  was  found 
to  converge  with  only  about  half  the  accuracy  achieved  with  the  upstream  differencing. 
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Table  2.  Solution  Procedure  for  the  Unsteady,  Inverse 
Boundary-Layer  Method 

Step  Procedure 

1.  Input  initial  values  of  H(x)  and  &*(x)  distributions  (flat  plate 
approximations)  and  ue(x)  distribution  from  inviscid  solution. 

2.  Evaluate  all  other  parameters  using  H.  Me,  and  Re0  at  every  spatial 
station. 

3.  Solve  equations  for  dd/dl  and  3ue/3t  at  every  spatial  station. 

4.  Compute  Me  and  H  from  values  of  ue,  0,  and  5*  using  backward 
differencing. 

5.  Use  Mc  and  H  distributions  as  approximations  for  next  time  step. 

6.  Go  to  step  2  and  do  until  convergence  (where  80/dt  and  3ue/dt  —  0). 

4.0  DIRECT  AND  INVERSE  UNSTEADY  BOUNDARY-LAYER  RESULTS 

To  choose  the  best  unsteady  method  for  coupling  with  the  Euler  code  to  obtain  an 
unsteady  viscous-inviscid  boundary-layer  code,  both  the  direct  and  inverse  unsteady 
methods  were  used  to  compute  attached  and  separated  boundary-layer  flows.  For  attached 
flows  the  inverse  code  was  run  in  the  pseudodirect  mode  to  use  the  same  experimental 
pressure  distributions  as  the  direct  method.  This  allows  a  more  direct  comparison  between 
the  two  methods.  For  separated  flows,  the  direct  method  uses  the  experimental  pressure 
distribution  while  the  inverse  method  uses  the  measured  displacement  thickness  distribution. 

The  two  attached  flow  cases  considered  (designated  Case  6  and  Case  9)  are  taken  from 
Ref.  17  which  includes  boundary-layer  data  for  an  RAE  2822  airfoil  at  several  Mach 
numbers  and  angles  of  attack.  Case  6  and  Case  9  upper  surface  measured  pressure 
distributions  were  used  as  input  for  both  codes  and  solutions  were  obtained  for  the 
boundary  layers.  The  airfoil  is  equally  divided  into  150  increments  in  the  present  study  and 
the  pressure  distribution  was  interpolated  for  each  point  to  begin  the  solution.  One  hundred 
fifty  increments  have  been  found  to  be  sufficient  for  describing  a  shock  located  on  the  air¬ 
foil.  Both  flow  cases  considered  include  a  shock  about  mid-chord.  Case  9  has  the  stronger 
shock.  The  codes  provide  a  continuous  solution  through  the  shock  region  and  no  turbulence 
enhancement  is  made  across  the  shock. 

The  various  boundary-layer  parameter  distributions  of  Case  6  for  both  the  direct  and 
inverse  solutions  are  given  in  Fig.  3.  The  two  solutions  are  indistinguishable  and  both  are 
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RAE  2822  Airfoil  (Upper  Surface),  Case  6 


M  =0.725 


Sym 


a  =  2,92  deg 
Re/c  =  6.5  x  10® 


O  Experiment 

-  Unsteady  Boundary- 

Layer  Code 

(Direct  and  Inverse 

Methods) 


b.  Skin  friction 

Figure  3.  Comparison  of  boundary-layer  parameter  measurements 
with  the  unsteady  bounday-layer  code  for  botb  direct  and 
inverse  methods  (Case  6). 
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plotted  as  one  line.  The  same  pattern  for  Case  9  is  shown  by  the  solutions  presented  in  Fig. 
4.  In  both  cases  the  shock  location  is  clearly  evident  as  an  abrupt  change  in  the  boundary- 
layer  parameters.  When  compared  with  available  boundary-layer  data  (Figs.  3  and  4),  the 
two  techniques  compare  well  upstream  of  the  shock  but  fall  below  the  data  downstream  of 
the  shock.  From  experience  with  the  steady  flow  codes,  this  result  is  not  unusual  when 
strong  shocks  are  present.  A  possible  reason  for  the  discrepancy  is  that  the  shock  may 
amplify  the  turbulence  of  the  flow  which  is  not  considered  in  either  the  shear  work  integral 
or  the  skin  friction  equation  of  the  present  method. 

RAE  2822  Airfoil  (Upper  Surface),  Case  6 

M  =  0.725  Sym 

CO  - - - 

a  =  2.92  deg  O  Experiment 

Re/c  =6.5  x  10®  -  Unsteady  Boundary- 

Layer  Code 


c.  Displacement  thickness 
Figure  3.  Concluded. 

Two  cases  of  low-speed  separated  flow  data,  compiled  by  Simpson,  et  al.  (Refs.  18  and 
19),  were  used  to  produce  solutions  by  both  the  direct  and  inverse  methods.  The 
experimental  displacement  thickness  distribution  was  used  as  input  to  the  inverse  method. 
Comparisons  of  both  solutions  with  data  are  presented  in  Figs.  5  and  6.  Flow  separation  (Cf 
=  0)  was  predicted  by  the  inverse  method  for  both  cases  of  Simpson  which  was  in  agreement 
with  the  measured  location.  The  direct  method  did  not  predict  flow  separation.  Solutions  by 
the  inverse  method  do  not  agree  well  with  momentum  thickness  data  downstream  of  the 
separation  point.  This  is  probably  because  the  shape  factor  correlations  used  are  applied 
beyond  their  useful  limits  (that  is,  H  >  4  occurs  in  the  present  case). 
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RAE  2822  Airfoil  (Upper  Surface),  Case  9 

M  =  0.734  Sym 

00  - 

ct  =  3.19  deg  O  Experiment 

Re/c  =  6.5  x  106  -  Unsteady  Boundary- 

Layer  Code 
(Direct  and  Inverse 
Methods ) 


b.  Skin  friction 

Figure  4.  Comparison  of  boundary  layer  parameter  measurements 
with  the  unsteady  boundary  layer  code  for  both  direct 
and  inverse  methods  (Case  9). 
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RAE  2822  Airfoil  (Upper  Surface) ,  Case  9 

Mm  =  0.734  Sym 

a  =  3.19  deg  O  Experiment 

Re/c  =  6.5  x  10®  -  Unsteady  Boundary- 


c.  Displacement  thickness 
Figure  4.  Concluded. 
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Sym  Simpson  I  Flow 

O  Experiment 


Unsteady  (Inverse  Method) 


a.  Momentum  thickness 

Figure  5.  Comparison  of  boundary-layer  parameter  measurements  of 
separated  flow  with  the  unsteady  boundary-layer  code  for 
both  direct  and  inverse  methods  (Simpson  1). 
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Sym  Simpson  II  Flow 

O  Experiment 

-  Unsteady  (Inverse  Method) 


x,  ft 

a.  Momentum  thickness 

Figure  6.  Comparison  of  boundary-layer  parameter  measurements  of  separated 
flow  with  the  unsteady  boundary-layer  code  for  both  direct  and  inverse 
methods  (Simpson  II). 
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Sym  Simpson  II  Flow 

O  Experiment 

-  Unsteady  (Inverse  Method) 

-  Unsteady  (Direct  Method) 


4  6  S  10 

x,  ft 

b.  Skin  friction 
Figure  6.  Concluded. 
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Two  important  points  are  made  from  the  unsteady  boundary-layer  results.  First,  both 
the  direct  and  inverse  method  produce  comparable  solutions  for  attached  flows.  Second,  the 
inverse  method  will  compute  a  separated  flow  solution,  whereas  the  direct  method  will  not. 
The  results  justify  coupling  the  inverse  boundary-layer  method  with  an  Euler  equation  code 
to  create  an  unsteady  interaction  code  capable  of  solving  flows  with  separation.  This  task  is 
the  focus  of  the  present  work. 

5.0  TRANSONIC  VISCOUS-IN  VISCID  INTERACTION  METHOD 
5.1  INTERACTION  PROCESS 

The  successful  application  of  the  inverse  unsteady  boundary-layer  code  in  conjunction 
with  an  Euler  code  would  greatly  enhance  the  progress  being  made  to  compute  three- 
dimensional  viscous-inviscid  interacting  flows.  The  applications  would  verify  that  the 
unsteady  interaction  methodology,  which  is  the  technique  to  be  used  in  three-dimensional 
computations,  is  feasible. 

The  interaction  process  used  in  the  steady  interaction  code  is  the  same  one  employed  by 
the  unsteady  formulation.  The  overall  interaction  scheme  is  based  upon  computing  the  Euler 
equations  for  several  time  iterations,  then  employing  the  viscous  method,  and  continuing 
this  viscous-inviscid  loop  until  steady  state  is  obtained.  The  unsteady  Euler  equations  are 
hyperbolic;  thus,  every  point  in  the  flow  field  influences  other  points  via  the  domain  of 
influence  in  the  time  plane  as  dictated  by  the  characteristics.  This  property  aids  in 
distributing  the  effects  of  the  boundary  layer  to  the  rest  of  the  flow  field. 

Whitfield’s  steady  interaction  method  applies  the  surface  source  model  of  Lighthill  (Ref. 
20)  as  a  coupling  device  between  the  viscous  and  inviscid  computational  regions.  This  choice 
permits  the  inviscid  computational  mesh  to  remain  fixed.  The  surface  source  model  imposes 
a  mass  flux  on  the  airfoil  of  the  form 

(ev)n  =  (eeM*),  (30) 

where  (ev)n  is  the  local  mass  flux  normal  to  the  surface  and  numerically  simulates  the 
displacement  effect  caused  by  the  boundary  layer. 

Physically,  the  surface  source  term  may  be  compared  to  a  blowing  effect  which  mimics 
the  boundary-layer  displacement.  There  continues  to  be  a  conservation  of  mass  within  the 
computational  domain  by  solution  of  the  continuity  equation.  Mass  entering  from  the  body 
is  a  boundary  influx  and  a  momentum  source,  as  at  the  entrance  plane,  and  is  matched  by  an 
additional  outflux  of  mass  at  the  flow  exit  plane.  The  added  mass  from  the  airfoil  displaces 
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the  effective  body  for  the  inviscid  flow  calculations.  However,  the  net  mass  added  in  the 
boundary  layer  and  wake  is  not  zero  and  a  global  mass  balance  is  not  satisfied.  Global 
momentum  considerations  then  show  that  a  thrust  must  be  exerted  on  the  airfoil.  Since  some 
of  the  added  mass  in  the  boundary  layer  is  removed  by  sinks  in  the  wake,  the  total  mass 
addition  is  small.  The  scheme  adopted  here  is  similar  to  the  method  in  which  the 
displacement  thickness  is  added  to  the  airfoil  for  each  new  inviscid  calculation.  However, 
the  present  method  has  the  advantage  that  the  computational  mesh  never  has  to  be  altered. 

After  the  flux  term,  (ev)„,  is  evaluated  in  each  pass  through  the  boundary-layer  routine, 
it  is  transferred  to  the  Euler  code  as  an  updated  boundary  condition  on  the  body  surface. 
During  solution  of  the  Euler  equations,  the  flux  term  is  held  constant  until  the  next 
boundary-layer  calculation.  Since  the  computations  were  made  for  airfoils,  the  upper  and 
lower  surfaces  were  computed  separately  within  the  boundary-layer  routine. 

In  each  viscous-inviscid  cycle  the  displacement  thickness  distribution  (initially  assumed 
as  a  flat  plate  approximation)  is  updated  by  information  from  the  previous  viscous-inviscid 
cycle.  This  is  done  after  the  Euler  equations  are  computed  for  several  time  cycles  and  before 
the  boundary  layer  is  computed  again.  The  method  used  to  specify  the  updated  5* 
distribution  is  Carter’s  method  (Ref.  21),  which  may  be  written  as 

6*(n+l>  =  5*<n)  +  aj(ue,v/ue,i-  1),  (31) 

where  £*(n  +  1)  is  the  updated  displacement  thickness;  5*(n)  is  the  displacement  thickness 
from  the  previous  viscous-inviscid  cycle;  uc>v  is  the  local  edge  velocity  obtained  from  the 
latest  boundary-layer  solution;  ue  i  is  the  velocity  obtained  from  the  previous  Euler  equation 
solution;  and  w  is  the  relaxation  parameter.  Details  of  the  boundary  conditions  in  the  Euler 
code  are  discussed  in  Section  5.2.  The  entire  solution  cycle  of  the  Euler  code  and  the 
boundary-layer  routine  is  continued  until  the  5*  distribution  is  converged.  The  term  u^/u^i 
which  approaches  unity  during  convergence  is  monitored  as  the  convergent  parameter.  An 
interaction  code  of  this  type  may  produce  better  solutions  than  an  integral  boundary-layer 
code  alone,  especially  if  shocks  are  present  in  the  flow.  While  the  flow,  including  the  shock, 
is  being  developed  by  the  Euler  code,  the  boundary  layer  is  continuously  interacting  and 
periodically  being  updated,  which  allows  the  boundary  layer  to  be  coupled  with  the  outer 
flow. 

The  unsteady  integral  boundary-layer  code  presented  in  Section  3.0  was  essentially 
exchanged  with  the  steady  boundary-layer  portion  of  the  steady  interaction  code,  creating 
an  unsteady  interaction  code.  After  the  exchange  was  accomplished  some  alterations  were 
needed  to  successfully  run  the  unsteady  interaction  code.  In  the  steady  interaction  code,  the 
incompressible  shape  factor  (H)  was  computed  directly  by  the  steady  boundary-layer 
solution  as  was  the  edge  Mach  number.  Since  all  correlations  are  based  upon  these  two 
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parameters,  any  other  parameter  required  in  the  solution  process  was  easily  obtained.  After 
each  boundary  layer  was  completed,  the  H  distribution  was  saved  and  used  in  the  next 
viscous-inviscid  cycle  as  a  new  approximation.  The  Mach  number  distribution  was  obtained 
from  the  Euler  solution  before  each  boundary-layer  calculation. 

i 

However,  the  parameters  computed  directly  by  the  unsteady  boundary-layer  solution  are 
0  and  edge  velocity  (ue).  After  each  viscous  cycle,  the  0  distribution  is  saved  for  the  next 
viscous-inviscid  interaction.  The  H  distribution  is  determined  from  the  computed  0  and  the 
updated  5*  distribution  (obtained  by  Carter’s  method).  After  the  H  distribution  has  been 
determined,  the  other  parameters  required  for  the  viscous  solution  process  may  be  computed 
from  the  shape  factor  correlations  and  edge  Mach  number  as  before. 

Another  problem  associated  with  the  unsteady  interaction  code  is  the  time-stepping 
iterations  of  the  viscous  part  of  the  code.  The  steady  viscous  solution  is  marched  spatially 
along  the  body,  while  in  the  unsteady  viscous  method,  a  number  of  time  iterations  must  be 
completed  within  each  boundary-layer  calculation.  The  most  economical  process  for  the 
present  cases  was  to  begin  with  a  few  iterations  (usually  50)  within  the  boundary-layer 
solution  and  arbitrarily  increase  the  number  of  viscous  iterations  with  each  viscous-inviscid 
loop.  The  last  time  the  boundary  layer  was  computed,  the  number  of  viscous  iterations  was 
arbitrarily  set  at  500  time  steps.  This  made  a  total  of  4,000  to  5,000  time  steps  within  the 
viscous  part  of  the  code.  Also,  a  total  of  approximately  1,000  time  steps  was  required  in  the 
inviscid  part  of  the  code.  If  no  boundary  layer  were  computed,  the  Euler  code  still  requires 
approximately  1,000  iterations  for  convergence  of  the  flow  cases  presented  here. 

5.2  EULER  EQUATION  CODE 

The  Euler  code  used  in  the  present  work  in  conjunction  with  the  inverse  unsteady 
boundary-layer  code  was  developed  by  Jameson  and  is  fully  described  in  Ref.  6.  The  code  is 
written  in  full  conservation  form  and  allows  the  representation  of  shocks  and  other 
discontinuities  in  flows  over  airfoils.  Far  field  boundary  conditions  allow  either  subsonic  or 
supersonic  flows.  All  solid  surface  boundary  pressures  are  extrapolated  from  the  flow  field 
values.  When  effects  of  the  boundary  layer  are  added  to  the  Euler  code,  the  flux  term 
associated  with  the  growth  in  the  displacement  thickness  is  included  at  the  solid  surface.  The 
method  of  solving  the  time-dependent  Euler  equations  is  a  four-equation,  fourth-order 
Runge-Kutta  scheme  with  local  time  stepping  and  central  spatial  differencing. 

A  C-mesh  is  used  in  the  Euler  code  (Fig.  7)  because  it  is  considered  best  for  airfoil 
solutions  in  which  the  wake  region  is  computed  because  it  allows  a  dense  spacing  in  the  wake 
region.  The  spacing  of  the  mesh  for  the  two-dimensional  cases  considered  in  the  present 
work  was  127  by  31  with  45  points  on  each  of  the  airfoil  upper  and  lower  surfaces.  The  mesh 
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Figure  7.  Representalive  sample  of  a  C-type  mesh  about  an  airfoil. 

may  also  be  aligned,  if  desired,  with  possible  discontinuities  such  as  shocks  or  slipstreams  in 
the  flow  Field.  No  changes  in  the  inviscid  mesh  are  required  for  applying  the  viscous  effects 
of  the  boundary-layer  code,  and  no  interpolations  are  needed.  The  mesh  remains  fixed  in 
time.  These  features  will  be  significant  factors  when  the  viscous-inviscid  method  is  extended 
to  three-  dimensional  problems. 

Viscous  solutions  of  the  wake  region  in  viscous-inviscid  iterations  are  computed  by 
establishing  a  small  value  (10-6)  for  the  skin  friction  and  applying  the  boundary-layer  code 
for  the  wake  as  an  extension  to  the  airfoil  surface.  The  wake  boundary  condition  is 
linearized;  i.e.,  is  applied  along  the  centerline  of  the  mesh  downstream  of  the  airfoil. 
Although  the  boundary  condition  is  not  applied  at  the  real  wake  location,  the  slipstream 
camber,  which  results  from  the  solution  of  the  code,  is  computed  at  the  correct  location.  A 
fair  representation  of  the  wake  may  be  made  in  this  manner. 

Boundary  conditions  required  in  the  Euler  code  include  specific  inflow  and  outflow 
conditions  for  either  subsonic  or  supersonic  flows  at  the  far  field  boundary.  Inflow 
conditions  for  supersonic  flow  are  set  to  free-stream  values  and  outflow  parameters  are 
extrapolated  from  the  interior  values.  For  subsonic  inflow  conditions,  the  density  is 
determined  by  extrapolation  while  other  properties  of  the  flow  are  evaluated  from  the 
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density  and  stagnation  conditions.  Pressure  is  set  to  the  free-stream  value  as  a  subsonic 
outflow  condition  while  velocities  are  obtained  by  extrapolation.  All  other  flow  parameters 
may  be  determined  when  these  far  field  properties  are  set. 

S3  RESULTS  OF  THE  STEADY  AND  UNSTEADY  VISCOUS-INVISCID  INTERACTION 

CODES 

Results  of  the  unsteady  interaction  code  developed  in  the  present  work  and  results  of  the 
steady  interaction  code  of  Ref.  4  are  discussed  and  compared  with  data.  Three  sets  of  data 
are  considered.  Case  6  and  Case  9  (Ref.  17)  of  Section  4.0  were  computed  with  both 
interaction  codes.  Another  data  case  from  the  same  reference  (denoted  Case  3)  is  also 
considered.  The  results  are  presented  in  two  parts.  The  first  includes  comparisons  of 
pressure  coefficient  data  with  computations  (Figs.  8,  9,  and  10);  the  second  contains 
comparisons  of  the  boundary-layer  data  with  computations  (Figs.  11,  12,  and  13).  Since  no 
lower  surface  boundary-layer  data  were  reported,  only  the  Case  9  lower  surface 
computations  are  presented  (Fig.  14)  to  show  that  the  trends  are  reasonable. 

For  all  three  cases  the  point  of  origin  of  the  boundary  layer  is  at  18-percent  chord.  For 
these  highly  loaded  aft-end  airfoils  it  is  essential  that  the  effective  camber  over  the  last 
20  percent  or  so  of  the  airfoil  be  estimated  correctly.  The  fixing  of  the  boundary  layer  at 
18-percent  chord  does  not  allow  the  boundary  layer  to  adjust  to  stagnation  point  movement. 
In  addition,  the  lower  surface  boundary  layer  may  not  be  predicted  too  well  in  the  strong 
adverse  pressure  gradient  downstream  of  50-percent  chord. 

The  predicted  pressure  distribution  for  Case  3  is  compared  to  data  for  both  the  steady 
and  unsteady  interaction  code  results  in  Fig.  8.  Both  computations  produce  essentially  the 
same  distributions  and  are  presented  as  a  single  curve.  The  comparison  with  data  is 
considered  good  except  on  the  upper  surface  ahead  of  the  10-percent  chord  position.  Such  a 
discrepancy  occurs  with  many  airfoil  calculations  and  is  related  to  the  lack  of  mesh 
resolution  near  the  leading  edge.  The  lower  surface  pressures  are  slightly  underpredicted.  If 

i 

the  flow  is  supercritical,  the  introduction  of  the  boundary  layer  at  18-percent  chord  causes  a 
jump  in  the  pressure  distribution  at  that  point  which  can  result  in  small  changes  in  lift 
attributable  to  the  differential  boundary-layer  growth  resembling  a  camber  change  to  the 
airfoil.  There  is  also  a  thickness  discontinuity  at  the  location  at  which  the  boundary  layer  is 
introduced. 
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Case  6  and  9  pressure  coefficient  comparisons  are  shown  in  Figs.  9  and  10.  The 
difference  between  the  steady  and  unsteady  calculations  is  apparent.  Both  methods  provide 
reasonably  good  results  when  compared  to  data,  and  only  a  slight  shift  in  the  shock  locations 
between  the  two  methods  is  evident.  Again,  the  largest  discrepancy  occurs  in  the  pressure 
distributions  of  the  lower  surface  and  the  initial  10-percent  of  the  upper  surface.  The  main 
point  to  be  observed  is  that  the  steady  and  unsteady  interaction  codes  agree  well.  In  all  three 
cases,  the  two  codes  give  equally  good  surface  pressure  predictions. 

Comparisons  of  the  computed  and  measured  upper  surface  boundary-layer  quantities 
for  Cases  3,  6,  and  9  are  presented  in  Figs.  11,  12,  and  13.  The  parameters  of  interest  are  the 
displacement  thickness  (6*),  momentum  thickness  (6),  and  local  skin  friction  (Cf).  The 
pressure  distribution  calculations  are  shown  as  an  aid  in  locating  the  shock  position.  Case  3 
and  6  computations  (Figs.  11  and  12)  show  good  agreement  between  both  methods  for  all 
three  parameters.  In  the  last  20-percent  chord  in  both  cases,  computations  and  data  have  a 
larger  discrepancy  than  at  the  upstream  positions.  Computations  of  Case  9  (Fig.  13)  provide 
the  best  agreement  with  the  data.  Near  the  trailing  edge,  good  predictions  are  made  by  both 
the  steady  and  unsteady  interaction  codes.  A  strong  shock  is  predicted  by  theory,  but  not 
enough  data  points  are  available  in  the  boundary  layer  to  determine  the  strength.  However, 
the  predicted  shock  location  and  strength  are  in  good  agreement,  Fig.  10.  As  shown  by  the 
skin  friction  curves,  both  computations  predict  a  near  separation  region  just  aft  of  the  shock 
and  near  the  trailing  edge.  The  two  computational  methods  agree  well  with  each  other. 

Several  assumptions  made  in  either  method  may  have  resulted  in  the  differences  denoted 
between  data  and  the  analytical  solutions.  Assuming  that  the  data  are  not  influenced  by 
extraneous  factors,  the  integral  method,  which  eliminates  turbulence  structure  in  the 
momentum  equation  and  the  simplistic  turbulence  model  used  in  the  shear  work  integral 
term,  may  be  contributing  factors  to  the  differences.  However,  the  overall  results  are  good. 
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RAE  2822  Airfoil,  Case  6 
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Unsteady  Interaction 
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Figure  9.  Comparison  of  the  measured  pressure  distribution  for  Case  6 
with  the  unsteady  and  steady  interaction  codes. 
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RAE  2822  Airfoil  (Upper  Surface),  Case  3 

Moo  =  0.60  Sym 

a  =  2.57  deg  O  Experiment 

Re/c  =  6.3  x  106  -  Unsteady  Interaction 


a.  Momentum  thickness 


b.  Skin  friction 

Figure  11.  Comparison  of  measured  upper  surface  boundary-layer  parameters  with 
calculations  of  the  unsteady  and  steady  interaction  codes  (Case  3). 
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RAE  2822  Airfoil  (Upper  Surface),  Case  3 
U  =0.60  Sym 
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Re/c  =  6.3  x  106  -  Unsteady  Interaction 
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c.  Displacement  thickness 
Figure  11.  Concluded. 
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RAE  2822  Airfoil  (Upper  Surface) ,  Case  9 

Moo  "  0.734  Sym 

a  =  3.19  deg  O  Experiment 

6 

Re/c  =  6.5  x  10  -  Unsteady  Interaction 

-  Steady  Interaction 


a.  Momentum  thickness 

Figure  13.  Comparison  of  measured  upper  surface  boundary-layer  parameters  with 
calculations  ui  the  unsteady  and  steady  interaction  codes  (Case  9). 
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RAE  2822  Airfoil  (Upper  Surface),  Case  9 
M  =0.734  Sym 

a  =  3.19  deg  O  Experiment 

Re/c  =  6.5  x  10®  - . Unsteady  Interaction 

—  —  —  Steady  Interaction 


b.  Skin  friction 
Figure  13.  Continued. 
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RAE  2822  Airfoil  (Upper  Surface),  Case  9 
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RAE  2822  Airfoil 
U  =  0.734 

oo 

a  =  3.19  deg 
Re/c  =  6.5  x  106 


(Lower  Surface) ,  Case  9 
Sym 

■  ■ '  Unsteady  Interaction 
-  Steady  Interaction 


c.  Displacement  thickness 
Figure  14.  Concluded. 


43 


AEDC-TR-84-7 


6.0  CONCLUSIONS  AND  RECOMMENDATIONS 

The  objective  of  this  work  was  to  produce  and  validate  an  unsteady  viscous-inviscid 
interaction  routine  based  upon  a  two-dimensional,  unsteady  integral  compressible 
boundary-layer  method  devised  by  Whitfield  and  the  Euler  equation  code  developed  by 
Jameson.  The  unsteady  interaction  code  was  compared  with  experimental  data  and  with  a 
steady  viscous-inviscid  interaction  code  which  utilizes  a  steady  integral  boundary-layer 
method.  The  results  indicate  that  the  unsteady  interaction  method  gives  satisfactory 
answers. 

The  major  advantage  of  the  present  work  is  that  the  method  is  computationally  faster 
than  a  Navier-Stokes  solution  while  providing  good  engineering  solutions  for  transonic 
airfoil  flows.  Another  advantage  of  the  method  is  that  the  inviscid  and  viscous  parts  of  the 
code  both  use  the  same  computational  mesh  spacing  in  the  streamwise  direction. 

The  major  disadvantage  of  the  method  is  the  use  of  a  simplistic  turbulence  model  within 
the  calculations.  Although  this  is  one  reason  for  the  speed  of  the  method,  the  turbulence 
modeling  at  a  shock  location,  for  example,  could  be  improved. 

To  improve  upon  the  present  method,  the  time-dependent  terms,  which  are  set  to  zero, 
should  be  better  approximated.  This  may  lead  to  an  accelerated  convergence  and  more 
accurate  solutions.  Further  testing  of  the  interaction  code  should  also  be  done  to  include 
wider  ranges  of  Mach  number  and  incidence.  Also,  cases  in  which  separation  is  known  to 
occur  should  be  computed  to  provide  information  about  the  method’s  ability  to  predict 
separation. 

Finally,  the  method  should  be  extended  to  three  dimensions  to  determine  the  difficulties 
that  may  arise.  Once  accomplished,  a  three-dimensional  interaction  code  can  be  developed 
to  compute  viscous-inviscid  solutions  of  realistic  configurations  in  transonic  flows. 
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APPENDIX  A 

DERIVATION  OF  EQUATIONS  SOLVED  BY  UNSTEADY  DIRECT  METHOD 

All  flow  parameters  are  evaluated  at  the  boundary-layer  edge.  The  subscript  “e”  has 
been  eliminated  and  is  understood. 

Expanding  Eq.  (17)  and  multiplying  by  u  will  yield 

=  jicl  _  A  (eu20)  _  _JH _ 

3t  at  2  gu  3x  dx  gu  3t  q  3t 


=  brfCf.M*,**) 


(A-l) 


Expanding  Eq.  (21)  and  multiplying  by  2u  will  yield 
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Using  the  definitions  6*  =  Hg.0*  and  8e  =  H#e  8,  expanding  Eq.  (A-l)  and  collecting  terms 
gives 

But,  H  and  Hpe  are  functions  of  H  and  Me  from  Whitfield’s  correlations.  Using  the  chain 
rule  and  assuming  dMe/dt  =  0  gives  the  final  form  for  Eq.  (A-4): 
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For  the  final  equations,  solve  for  M/dt  and  3H/3t  between  Eqs.  (A-3)  and  (A-5)  to  give 
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APPENDIX  B 

DERIVATION  OF  EQUATIONS  SOLVED  BY  UNSTEADY  INVERSE  METHOD 
In  a  similar  fashion  as  in  Appendix  A,  expanding  Eq.  (17)  and  multiplying  by  u /0  results 


m 
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udt  0  dx  dx  2  6  6  d t  * 
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q  dx  q  dt  Cz 

Expanding  Eq.  (21)  and  multiplying  by  2u/0  results  in 


_L  39 

6  d t  + 


T'-> 

K  ♦  2  (*  -  4 )]  -S'  + 


89 

dx 


(B-I) 


-  u 


dx 


(B-2) 


_  _L  J_  8q  _  (1  +  Hy  -  Hgg)  se  _ 

9  8t  q  dx  e  dt 

This  pair  of  equations  may  now  be  solved  for  86/ dt  and  du/dt  by  Cramer’s  rule: 


Defining 
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then 


dfl  =  A22  C|  -  A12  C2  _  Aa  C\  -  A12  c2 

dt  An  A22  -  A]2  A2j  Aji  A22 


(B-3) 


du  _  An  c2  -  A2i  C]  _  c2 
3t  An  A22  -  Aj2  A2j  A22 


(B-4) 


These  equations  can  be  solved  as  long  as  the  denominator  (Ha./u0)  does  not  approach 
zero  (i.e.,  where  Ha*  —  O'  or  u0  —  00).  This  condition  will  not  apply  since  Ha*  *  0  and  ufl  is 
finite. 


Both  bi  and  b2  in  the  direct  method,  and  ci  and  c2  in  the  inverse  method  include 
derivatives  in  time  of  flow  parameters  other  than  the  ones  sought  in  the  solution  process.  It 
should  be  noted  that  the  time-dependent  terms  are  set  to  zero  in  the  formulation  to  solve  the 
equations.  The  end  result  is  to  obtain  a  steady-state  solution  in  which  these  terms  have  no 
influence.  However,  during  the  solution  process,  an  assumption  must  be  made  for  these 
terms. 
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NOMENCLATURE 

All,Ai2, 

A21 » A22  Denotes  elements  of  the  solution  matrix  (Appendix  B) 

bi,b2  Functions  defined  in  Appendix  A 

Cp  Pressure  coefficient,  (p-p»)/(l/2  q*) 

c  Airfoil  chord 

C|,C2  Functions  defined  in  Appendix  B 

Cf  Local  skin  friction  coefficient,  2r0/(eel$ 

Dg  Shear  work  integral  term,  defined  by  Eq.  (20) 

e  Energy,  defined  by  Eq.  (5) 

Fjj  Parameter  used  in  Navier-Stokes  equations,  defined  by  Eq.  (3) 

fy  Parameter  used  in  Navier-Stokes  equations,  defined  by  Eq.  (3) 

G  Parameter  used  in  Navier-Stokes  equations,  defined  by  Eq.  (2) 

H  Incompressible  shape  factor,  5V0 

H{»  Shape  factor  based  on  £*,  defined  by  Eq.  (24) 

H4*u  Shape  factor  based  on  £*,  defined  by  Eq.  (26) 

Hj*  Shape  factor  based  on  0*,  defined  by  Eq.  (25) 

Hje  Shape  factor  based  on  dQ,  defined  by  H#p  =  0e/0 

M  Mach  number 

p  Static  pressure 

q;  The  heat  flux  component 
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Qoo  Free-stream  dynamic  pressure,  gooU^/2 

Re/c  Reynolds  number  based  upon  airfoil  chord  length,  CooUooC/^o, 

Ree  Reynolds  number  based  upon  the  momentum  thickness,  Qeu.^/i^ 

t  Time 

u  Velocity  in  the  x-direction 

v  Velocity  in  the  y-direction 

x  Coordinate  along  body  surface 

xi>x2ix3  Arbitrary  orthogonal  coordinate  system 

y  Coordinate  normal  to  the  body  surface 

Greek  Symbols 

or  Angle  of  attack 

7  Ratio  of  specific  heats 

Kronecker  delta 

6*  Boundary-layer  displacement  thickness,  defined  by  Eq.  (13) 

5*  Boundary-layer  velocity  thickness,  defined  by  Eq.  (18) 

7j  Bulk  viscosity 

6  Boundary-layer  momentum  thickness,  defined  by  Eq.  (14) 

0*  Boundary-layer  energy  thickness,  defined  by  Eq.  (19) 

0e  Boundary-layer  density  thickness,  defined  by  Eq.  (15) 

X  Second  coefficient  of  viscosity 
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V 

G 

-<gu'v'> 

<6V)n 

<*j 

T 

ra 

<i> 

Subscripts 

e 

o 

00 


Molecular  or  first  coefficient  of  viscosity 
Kinematic  viscosity  coefficient 
Density 

Reynolds  stress 

Local  normal  mass  flux,  defined  by  Eq.  (30) 

Deviatoric  stress  tensor 

Total  shear  stress 

Denotes  the  viscous  stress  tensor 

Relaxation  parameter  in  Eq.  (31) 

Boundary-layer  edge  conditions 
Wall  conditions 
Free-stream  conditions 


Superscripts 

( — )  Low  speed  or  incompressible  value 
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